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Abstract 

Monitoring the health of ancient artworks requires adequate prudence because of the sensitive nature 
of these materials. Classical techniques for identifying the development of faults rely on acoustic testing. 
These techniques, being invasive, may result in causing permanent damage to the material, especially if 
the material is inspected periodically. Non destructive testing has been carried out for different materials 
since long. In this regard, non-invasive systems were developed based on infrared thermometry principle 
to identify the faults in artworks. The test artwork is heated and the thermal response of the different layers 
is captured with the help of a thermal infrared camera. However, prolonged heating risks overheating 
and thus causing damage to artworks and an alternate approach is to use pseudo-random binary sequence 
excitations. The faults in the artwork, though, cannot be detected on the captured images, especially if 
their strength is weak. The weaker faults are either masked by the stronger ones, by the pictorial layer 
of the artwork or by the non-uniform heating. This work addresses the detection and localization of the 
faults through a wavelet based subspace decomposition scheme. The proposed scheme, on one hand, 
allows to remove the background while, on the other hand, removes the undesired high frequency noise. 
It is shown that the detection parameter is proportional to the diameter and the depth of the fault. A 
criterion is proposed to select the optimal wavelet basis along with suitable level selection for wavelet 
decomposition and reconstruction. The proposed approach is tested on a laboratory developed test sample 
with known fault locations and dimensions as well as real artworks. A comparison with a previously 
reported method demonstrates the efficacy of the proposed approach for fault detection in artworks. 

Index Terms 
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Non-destructive testing, subspace decomposition, wavelet transform, fault detection, mutual informa¬ 
tion, Infrared thermography, PRBS excitation 

I. Introduction 

The detection and identification of faults in different materials has seen major emphasis over the 
last decades, both in research context and in industrial applications. Diverse techniques like Gamma 
and X-rays, ultrasounds, Foucault currents and Nuclear magnetic resonance are conventionally utilized 
for non-destructive testing of materials HI, IH, IS, IH, ||5l, IH, Q. Recently, techniques focusing on 
thermal radiations, such as photo-reflectance, photo-acoustics, mirage effect and photo thermal radiometry, 
have been tested on the laboratory scale for non-invasive testing of materials (H. One such technique 
is the photo-thermal radiometry (or photo-thermometry) which requires a relatively simple and low-cost 
experimental setup Q- The test material is excited by a flow of heat, resulting in a change in local thermal 
conditions of the material. The thermal response of the material to this excitation is then captured by a 
thermal infrared camera. The acquired thermal response depends on different parameters of the material 
such as thermal conductivity, diffusivity, emissivity and specific heat as well as the excitation used at the 
input. More specifically, the above properties manifest themselves in the thermal response depending upon 
different factors such as fissures or holes in the material, material structure, physicochemical processes 
taking place in the material, delamination, sedimentation, etc. The detection of faults in ancient artworks 
is a challenging problem requiring adequate precaution to identify the faults without deteriorating the 
material, even slightly. The conventional method used for exciting the material is the pulse excitation 
method which consists of heating the material over a relatively long time and then allowing it to cool 
down ifTOll . In the current research endeavor the frescos(Italian murals) are the subject of interest which 
are layered structures. Due to the environmental effects air gaps occur at the interface of the layers. This 
is the working definition of a defect in the artwork which will be used through out the current paper. The 
acquired thermal response is relatively easier to process from fault identification purposes. However, the 
prolonged heating involved in the process risks overheating and thus causing damage to the material. An 
alternate approach is to use random excitations like the Gaussian and pseudo-random binary sequence 
(PRBS) excitations. The step heating results in a greater transfer of thermal energy to the material under 
analysis as compared to the PRBS excitation.In photo-thermal radiometry the detection of fault is based 
on the variations in the thermal radiations based on the presence of faults hence a greater transfer of 
thermal energy will yield better contrast. While these techniques are more favorable from material health 
conservation perspective, they necessitate more intensive processing to retrieve useful information about 
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the defects. 

The useful information should reveal the localization maps of different defects in the material after 
having suppressed the effects of pictorial layer, measurement noise, potential interference patterns due 
to multiple excitation sources and inhomogeneous illumination of the test material. Recently, some 
signal processing based methods were deployed for detecting defects in artworks while using PRBS 
excitations ifTTI . Although, the method revealed some interesting results regarding the material defects, 
the interpretation of the final results was not very evident. In this paper, we propose a new processing 
algorithm based on wavelet decomposition for identification of defects. The algorithm exploits the 
multi-resolution capability of wavelet transform and approximates an efficient subspace, which contains 
information on the defects, while removing other influences from the acquired thermal response. The 
proposed subspace decomposition approach is applied spatially and the resulting subspaces for different 
acquisitions are then temporally processed to extract a single representative subspace for the defects. Since 
there is no universal criterion for the selection of wavelet basis function, we propose a selection criterion 
as well as an automated mechanism for determining the wavelet decomposition level. Even though the 
main focus of the paper is proposition of a defect detection algorithm, it is demonstrated that the intensity 
of the final detection parameter can be associated to the intensity of the defect. The proposed method 
is validated on a sample test material, containing defects of various depths and diameters. The detection 
parameter allows localizing these faults in the test sample with different intensities. The method is also 
applied on a real artwork to reveal interesting information about the defects. The paper is organized to 
start with a presentation of the experimental setup, the test samples and the acquired data set. This is 
followed by the proposed methodology section comprising the mathematical formulation of the problem, 
the details of the proposed algorithm and the different parameters involved therein. The results of the 
proposed approach on a test sample are discussed in the subsequent section. The paper is finally concluded 
fo highlighf fhe pros and cons of fhe proposed approach along wifh an insight into some future work in 
the domain. 


II. Experimental Setup and Raw Data Set 

In this section, the experimental setup for photo-radiometry system including the excitation and acqui¬ 
sition process is described, along with the presentation of a sample data set. 
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(a) Experimental setup demonstrating sample illumina- (b) A typical pseudo-random excitation sequence 

tion and recording 


Fig. 1. Experimental setup for the non-invasive testing of a sample material. The material is excited by two 500 W halogen 
lamps, simultaneously controlled with a pseudo-random excitation sequence. The resulting thermal response of the material is 
captured by a thermal long wave infrared camera. 


A. Excitation and Acquisition System 

The experimental setup for photo thermometry is illustrated in Figure |l(a)[ The artwork is heated 
simultaneously hy two halogen lamps, each one with a power of 500 W. The lamps are symmetrically 
placed around the infrared camera and placed at approximately 50 cm from the artwork. The temporal 
excitation pattern of these two lamps is controlled hy a pseudo random binary sequence (PRBS). The 
smallest duration of each binary pulse in this sequence represents the excitation time (Tg). The experiment 
described in this paper is carried out for different excitation times as represented in Table JlThis was 
done to demonstrate the generality of the proposed algorithm across different experimental parameters. 
An illustrative excitation pattern is shown in Figure |l(b)[ 

The incident flow of heat strikes the surface of the material and its subsequent absorption depends on 
the physical properties of the material. The surface defects, such as cracks, might appear instantly due to 
the different optical properties and require other methods for detection. The sub-surface defects appear 
after some time depending on their respective depths. The major challenge is posed by the defects that 
are not visible at the surface and are located at different or varying depths. Such defects cause air pockets 
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TABLE I 

Excitation TiMES(Te) and the corresponding number of pulses in the PRBS sequences ( Nt ) 


Te 

0.5s 

Is 

2s 

5s 

10s 

Nt 

2048 

1024 

512 

256 

128 


in otherwise uniform material, resulting in different distributions of the heat inside the materials. The 
radiated thermal energy is modified and presents an opportunity to identify the faults and to quantify 
them. In brief, different defects into a material exhibit different thermal responses to the same excitation. 

The thermal response of these different areas subjected to the excitation sequence are then captured 
by a A20 FLIR thermal Infrared (IR) camera. This bolometer camera, working in long wavelength range 
(7.5 to 13 //m) and having a thermal sensitivity of 0.12°C' at SO”^, was placed around 50 cm from the 
sample. The excitation of the lamp and the IR camera were synchronized, the former being in a slave 
mode. A thermal image is obtained for each pulse duration of the input PRBS sequence, the total number 
of such images is thus same as the number of pulses in the PRBS sequence, Nf. The number of pulses 
in the PRBS sequences corresponding to the different excitation times is mentioned in Table U 

Even though, the relative positions of the halogen lamps and the camera have been carefully setup 
after exhaustive testing, the problem of inhomogeneous illumination can not be ruled out. It will be 
demonstrated that the proposed wavelet based subspace decomposition approach allows estimation of a 
background affect which may be attributed to inhomogeneous illumination. The preceding description of 
the experimental setup indicates that the thermal behavior of the material is spatially non-stationary. The 
subsequent sections will reveal how the proposed approach allows exploiting this spatial non-stationarity. 

B. Test Materials 

Before developing a mathematical formulation for the problem of defect detection, we briefly shift our 
focus to the presentation of the test samples used in this work. In particular, we consider two types of 
test materials, one with faults of different dimensions simulated by introducing holes in a plaster material 
and other are real artworks containing faults of varying geometry and dimensions. 

1) Laboratory test sample : The complete acquisition chain of the experimental setup will be demon¬ 
strated with the help of this test material to enhance understanding of the underlying problem, which 
needs to be modeled. It will also serve as a ground truth against which different parameters of the 
proposed algorithm will be established. The test material is composed of a white plaster substrate with 
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a thickness of approximately 20mm (Fig. |2l). The front side of the material, which is exposed to the 
excitation, does not contain any significant defects, hut still has some minor irregularities, especially 


toward the edges (Fig. 2(a) i. As descrihed earlier a defect in an artwork is typically characterized hy the 
slight delamination of the material in the vicinity of the defect causing air pockets within the material. 
The hack side of the material was therefore modified fo create defects of different dimensions (Fig. |2(h)| ). 
In particular, differentiation amongst these defects arise from differences in their depths and diameters. 
Twelve different holes were created in the material and their dimensions (depths and diameters) are 
marked in Fig. |2(h)| and also summarized in Table HIl It is pertinent to mention that these depths are 


measured from the front side of the material. Thus, hole 1 (Fig. |2(h)| top right) is at a depth of 4 mm 
from the front surface. The holes have been numbered in sequence starting from top right (hole 1) and 
moving left and down to bottom left (hole 12). As illustrated in Fig. |2(b)[ three different depths, (4 
mm, 6 mm, and 8 mm) with four different diameters (3 mm, 6 mm, 8 mm, and 10 mm) for each of 
these depths, have been tested in this experiment. It is assumed that each fault has uniform depth and 
diameter signifying that at the spatial location of the fault, there’s an air pocket across the material up 
till the depth of the given fault. As the surface is heated the thermal energy gradually diffuses through 
the depth of the material. The rate of the diffusion depends on the thermal diffusivity of the material. 
The thermal diffusivity of the plaster of Paris (the main constituent of the fresco) is greater than that of 
the air (present at the fault locations). This results in a difference in the rate of diffusion of heat and in 
result the amount of back scattered radiation results in detection of fault regions(as they present higher 
thermal intensity). Shallower and wider faults offer bigger air pockets which effect the thermal diffusivity 
and conductivity. This translates physically into the 4 mm hole (depth) more likely to depict a higher 
heat intensity as compared to the 6 mm and 8 mm holes. The overall impact on the outcome is that the 
acquired thermal response is more intense at the location of shallower and wider holes as compared to 
deeper and narrower holes. 


Photo-radiometry is an inspection technique, which works at depths closer to the contact surface 0- 
While, the higher intensity faults may be identifiable in the raw thermal images, the smaller faults are 
often masked. The goal of this work is to extract such faults from the background. The test sample 
described above will also allow studying the influence of the diameter (spread) of a certain defect for a 
given depth. Moreover, it will reveal the impact of processing raw data through the proposed algorithm 
in terms of enhancing the detectability of the faults at different depths. It should be mentioned that real 
defects in an artwork will not necessarily be circular in shape, but the purpose of the current research is 
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(a) Front side directly exposed to excitation (b) Back side containing holes to simulate defects 


Fig. 2. Front and back sides of the test material prepared in laboratory. The material is constructed from plaster with the front 
side exposed to the flow of heat containing no holes or other significant anomalies. The back side of the material consists of 
12 holes of varying depths and diameter to simulate defects. Depths are indicated from the front side. 


TABLE II 

Characterization of laboratory prepared test sample (Fig.IH): co-ordinate mapping of the faults’ 
(HOLES’) centres; DIMENSION OF THE FAULTS; AND THE IMPROVEMENT IN THE SNR (DIFFERENCE BETWEEN THE SNR 
AFTER THE APPLICATION OF PROPOSED ALGORITHM AND THE SNR OF THE RAW DATA) AFTER THE APPLICATION OF 

PROPOSED ALGORITHM CALCULATED AS IN EQUATION ([Hi. 


Hole 

Coordinates(pixels) 

Diameter 

Depth 

Improvement in SNR via proposed algorithm 


iVx 

Ny 

(mm) 

(mm) 

(dB) 






Te = 0.5s 

Te = Is 

Te = 2s 

Te = 5s 

Te = 10s 

1 

39 

167 

10 

4 

9.9579 

8.0364 

7.5142 

9.7085 

9.8548 

2 

39 

127 

8 

4 

4.5801 

3.8630 

5.2932 

9.2458 

3.5419 

3 

42 

83 

6 

4 

-1.5678 

-2.8933 

2.4168 

5.5109 

2.2004 

4 

41 

30 

3 

4 

9.3600 

7.9558 

7.3496 

-1.0940 

10.7252 

5 

77 

167 

10 

6 

4.9934 

4.5480 

3.7430 

6.2764 

6.1833 

6 

76 

127 

8 

6 

2.2338 

-2.6496 

2.8383 

5.4974 

6.4358 

7 

83 

83 

6 

6 

4.8059 

5.3750 

5.4407 

2.1080 

6.2786 

8 

86 

31 

3 

6 

-5.5033 

2.9049 

2.5581 

1.0447 

-0.1188 

9 

117 

169 

10 

8 

-5.1488 

-2.4098 

-2.8700 

2.4369 

1.4146 

10 

118 

127 

8 

8 

-2.1426 

1.6670 

-2.9004 

5.2986 

-5.9618 

11 

122 

83 

6 

8 

2.6452 

-2.6246 

0.4480 

-2.2460 

-1.4110 

12 

130 

31 

3 

8 

3.8080 

-0.5390 

-4.0005 

-2.5897 

5.8405 
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not to indulge in particular defect geometries but rather to establish the detectability limits. Nevertheless, 
the viability of the proposed approach for irregular defects will be demonstrated on real artworks in the 
results section. 

2) Real Artwork: In order to validate the proposed approach, thermal responses of a real artwork, 
namely Mural 1, were acquired. This artwork represents a replica of the “Saint Christopher with the 
Christ Child”, a Florentine fresco realized in the end of the 14*^ century, created in laboratory according 
to the technique of the Italian primitives ifTTl . Fig. |3(a)| shows a picture of Mural 1, along with the 
designated faults in this artwork (Fig. |3(b)[ ). Five inclusions of plastazote have been introduced in the 
fresco in its manufacturing process. These five faults are labeled A to E: A is tilted with varying depths 
of 3 to 10 mm and a thickness of 5 mm; B is at a depth of 5 mm with a thickness varying between 3 
mm and 10 mm; D is located at a depth of 3 mm and has a thickness of 3 mm; E is at a depth of 10 
mm, lying beneath D, and has a thickness of 5 mm; finally, C is a powdery defecl located af a depfh of 
3 mm and has fhickness of 5 mm. These faulls are of non-uniform shapes. 

Anofher real artwork used during fhe course of fhis work is named Mural 2. No information regarding 
fhe posifion of fhe sub-surface faulfs is known prior fo fhe application of fhe proposed algorifhm. The 
surface arfifacf is a piece of gold foil which in fhe course of presenf invesfigalion is nol of inferesf and 
will hence be removed using posf processing techniques as discussed later on. 

III. Proposed methodology: selective sub-space extraction eor defect detection 

In fhis secfion, we will presenf fhe proposed approach for fhe defection of defecfs. Before delving 
info fhe defails of fhe algorifhm, we will associafe a mafhemafical formulalion fo fhe acquired fhermal 
responses. Moreover, we will also briefly discuss a mefhod previously developed by fhe aufhors in order 
fo subsfanfiafe fhe merifs of fhe approach proposed in fhis paper. The acquired raw dafa consisf of Nt 
fhermal images, as shown in Fig. |4l Note fhaf images acquired on fhe fronf side were flipped on fhe y axis 
in order fo have fhe same positioning of fhe holes shown in fhe back side (see Fig. |2(b)| ). The raw dafa 
cube, y, exisf in fhe Hilbert space y € where and Ny correspond fo fhe dimension of 

fhe lesl material (pixels) as capfured by fhe camera and Nt corresponds fo fhe fofal number of acquisitions 
in time (samples). While and Ny will be fixed for a given fesf sample (assuming fhe experimenfal 
sefup is nol changed and fhe resolution of fhe camera is nol exceeded), fhe value of Nt will depend on 
fhe excilafion time (Tg) as represented in fhe Table Jl The raw dafa confains information from mulliple 
sources including defecfs, background or picforial layer, inhomogeneous illumination and measuremenl 
noise. The Hilbert space conlaining fhe raw dafa can be decomposed info ils subspaces, each represenling 
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(a) Picture of Mural 1 


(b) Fault map of Mural 1 


Fig. 3. A real artwork, Mural 1: (a) photograph of Mural 1 displaying no visible defects; (b) fault map highlighting the locations 
of faults and their irregular geometry for performance analysis. The faults A to E are incorporated with different dimensions. 


single or multiple afore mentioned sources. Subspace decomposition is a processing technique which 
allows more meaningful projections. In our case, we are looking to decompose the data into background, 
defect (useful) and noise subspaces. Subspace decomposition based matrix filtering techniques have been 
used in diverse applications l[T2ll . ifT^ . ifldll . ifTSl . ifT^ . Singular Value Decomposition (SVD) and Wavelet 
decomposition are common methods of sub-space decomposition. In this paper, the SVD based technique 
for defect identification will be compared with the proposed wavelet based decomposition technique. 

A. State-of-the-art decomposition of thermal responses 

An SVD based defect detection scheme in artworks using the thermal responses was developed by 
Vrabie et al. m. The scheme was based on the hypothesis that the evolution of thermal response over 
time is a distinguishing factor between defects and background. It involved rearranging the data cube, 
y G sfjWxiVa,xArt ^ 2-D data set Y € where Ng = N^xNy, represents the spatial dimension. 

This 2-D space was then decomposed into three subspaces, the first one representing background, the 
second one useful information containing defects and the third one the noise components. The number of 
singular values to be considered for each of the three subspaces was determined empirically. Specifically, 
only one singular value confaining 96% of energy was used for exfracting the background, whereas the 
so called useful subspace containing defects was synthesized from few subsequent components and the 
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Fig. 4. Raw data y acquired on the laboratory generated test material of Fig. One thermal image per acquisition pulse is 
acquired thus generating a data cube. The colors represent the thermal intensities as measured by the camera. 


lowest energy singular components were attributed to the noise. The second subspace was selected for 
further processing to detect the locations of defects. The higher order statistics (HOS) of skewness and 
kurtosis were then computed for the time series in the useful subspace to obtain the results in space, 
which were then rearranged to obtain an image map in for both kurtosis and skewness. It was 

demonstrated that these skewness and kurtosis maps can help in identification of defects. However, the 
interpretation of the results obtained with this SVD and HOS based approach was not systematic. 

Although, SVD is a matrix filtering technique, it was utilized in the temporal domain without any 
focus on explicitly exploiting the spatial patterns of the defects. The high energy first component of 
SVD contains information on the average behavior of the input data. Since, the SVD-HOS approach, 
discussed above, discards the first subspace, it therefore runs the risks of suppressing high energy defects 
along with the background. This fact is attributable to the poor separability of different components by a 
second-order statistics based approach such as SVD. The limitations of this SVD-HOS based approach 
led us to the development of an alternate approach for detection of defects. This multi-resolution analysis 
based approach also performs the subspace decomposition but in spatial rather than the temporal domain. 
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B. Wavelet based subspace decomposition 

The raw thermal response acquired hy the camera at a given instant, Yt G incorporates 

information on the defects hut this information is often masked hy many undesirable effects such as 
inhomogeneous illumination, background pictorial layer and measurement noise. The inhomogeneous 
illumination is caused by irradiation from two different sources whose placement, though carefully 
adjusted, could not make the captured image illumination invariant. The characteristics of the pictorial 
layer of the test material such as pattern and/or the paint also manifest themselves strongly in the acquired 
raw image. In addition, the defects can have different dimensions in terms of the wear-ability of the 
material. The stronger defects will have more significant intensities and will appear sooner as compared 
to the weaker ones because of the fact that for stronger defects (e.g. those closer to the surface), the 
heat diffusion will be greater as compared to weaker ones. The spatial pattern of the raw thermal image 
is therefore non-stationary. In order to better exploit this non-stationarity, we resort to the joint-time 
frequency analysis of wavelet transform. The underlying principle is that of the multi-resolution analysis 
whereby it is possible to decompose the input data into multiple subbands. 

Discrete Wavelet Transform (DWT) is a powerful tool for performing multi-resolution analysis and has 
been exploited in numerous fault detection applications |[T5l . ifTTl . lITSl . lfT9ll . EOll . lfT3l . At any level of 
decomposition ‘L’, DWT decomposes the input data into a high pass or detail subspace Dl and a low 
pass or approximation subspace A^, with input to the level L = 1 being Yt and to any subsequent level 
L, the approximation of the previous level Ai_i. To apply the DWT on 2-D datasets, the filtering is first 
applied along one dimension and then the resulting images are processed along the other dimension. Dl 
represents the detail components resulting from high pass filtering in at least one dimension (horizontal 
or vertical) whereas Al corresponds to the low pass filtering output in both dimensions. As depicted in 
Fig. [5l application of DWT on images yields 4 sub-images, out of which the sub-image resulting from 
the two low-pass filters represent the approximation subspace and the remaining images are summed up 
to get the detail subspace. The output of the wavelet decomposition process depends on the basis function 
used for computation. In general, the basis function should be similar to the features being searched for 
in the data, however, the exact signatures of the features being not known, this choice is mostly empirical. 
In this paper, we perform an exhaustive search approach to select the best basis for our defect detection 
application which will be discussed later. Another important parameter is the number of levels to be 
used for decomposition which will mainly depend on the separability between the strongest defects, the 
background pictorial layer and the illumination effects. We utilize a criterion based on regional mutual 
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Fig. 5. 2 — D Wavelet decomposition at level L. The low pass image Al-i at a given time instant is first horizontally filtered 

through 1-D low pass and high pass filters followed by a vertical pass on the results. Dl corresponds to the three details and 
Al represents the approximation space. 


information ll^ for the selection of number of decomposition levels. 

1) Proposed algorithm: The acquired data at a given instant, Yt, being a mixture of different sources 
(inhomogeneous illumination, measurement noise, pictorial layer, background effects and useful defects), 
the goal is to separate out the useful band pass region containing the information of the defects. DWT 
with its multi-resolution capability acts as a filter bank with pass bands of different widths. Since our 
objective is the localization of fault in the spatial plane, we apply the 2D wavelet transform in the spatial 
domain at each time instant. The proposed algorithm is presented in Fig. 0 and involves three principal 
stages: multi-level decomposition of each image, selective reconstruction and temporal averaging to obtain 
the final detection result in terms of a fault map (image). The relevance of the DWT output depends on 
the choice of the basis function. The signal representation using a basis set essentially corresponds to 
finding the inner product so to obtain the best representation, the basis function must be similar to the 
analyzed signal. The acquired data, Yt, is dependent on the excitation sequence, the thermal properties 
of the test material and the nature of existing faults. As the nature of the faults is not known beforehand 
and we are not physically quantifying the material, we propose a method for determining the best basis 
function for our application based on an exhaustive search approach. 

After decomposition using selected basis, the useful information containing the defects is extracted by 
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Fig. 6. The proposed scheme based on wavelet decomposition and selective reconstruction. The main parameters of the scheme 
are the wavelet basis (j) and the number of decomposition levels, L. Yr is reconstructed from intermediate level details and 
contains information about the faults. 


reconstruction using appropriate subspaces. The space containing the raw data Yt is decomposed into 
three subspaces as: 

Yt = Yb + Yr + Yd, (1) 

where Yb represents the subspace corresponding to the background pictorial layer and uneven illumination 
effects as obtained in the approximation Al at level L, Yr corresponds to the useful subspace containing 
fault information and is reconstructed from intermediate level details and Yd represents the discarded 
subspace corresponding to the high frequency noise appearing in the lower level details. 

Three important parameters in the proposed algorithm (Fig. O are the wavelet basis selection; the 
number of decomposition levels and; the selection of the details to reconstruct Yr. Before focusing on 
the selection criteria for these three parameters, we consider the third module of the overall processing 
algorithm (Fig. O. The first two modules of image decomposition and reconstruction are applied to every 
temporal acquisition thus resulting in a representative fault subspace at every time instant. Since the final 
goal is to localize the faults in the spatial domain while using PRBS excitations, the fault subspaces Yr, 
at different time instants, are temporally averaged to obtain the final result Ydet G Ydet thus 

represents the spatial fault map of the sample being tested. 

2) Parameter selection criteria: An important aspect in the application of wavelet transform is the 
choice of basis function to be used in the decomposition and reconstruction of the subspaces. Although, 
wavelet analysis has found widespread applications in diverse domains, it is not a data driven technique 
and there does not exist a global rule to select the best basis. The most frequent approach is to perform 
exhaustive testing with different basis functions and select the one which produces the best desired 
results. As the current problem is based on the classification of thermal image pixels into fault and non- 
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fault regions we require a ground truth to establish the effect of the values of different parameters on the 
classification accuracy of the proposed algorithm. The test sample created in the laboratory represented in 
the Figure |2(b)[ presents an excellent test subject as it simulates faults of differnt depths and diameters. 
For our application, we propose a criterion for best basis selection by formulating and subsequently 
minimizing an error function. We exploit the detection results on the test sample (Fig. |2ll to select the 
basis which provides best separation of the desired faults from the background without compromising 
on the weaker faults. The cost function is formulated as: 

E l|Ydet(win)i • - Ydet(win), ir 

_ ( 2 ) 

||Ydet(win)^^- - ^avg,j\\ 

where, $ represents the set of all the basis functions (j)j, which are tested, Ydet(win)-^ represents the 
detection parameter at the output of the algorithm computed in a window around the fault for the 
basis, ^avg,j represents the average value of the detection parameter in the background, where there are 
no known faults for the basis and ||.|| represents the matrix norm. The numerator of Eq. |2] ensures 
that the difference between the detection results for all the faults is minimized with respect to the most 
significant fault that corresponds to Hole 1, whereas the denominator maximizes the differences between 
the most significant Hole 1 and the background that corresponds to the noise. In fact, viewing the problem 
in terms of the classification of raw image into defective and non-defective regions, the minimization of 
this cost function ensures that the inter-class distance is maximized while simultaneously minimizing the 
intra-class distance. The best basis cpj is the one that provides minimum value for this cost function. 

It is pertinent to mention that the proposed approach has been formulated given a ground truth about 
the nature (depth and diameter) and the location of the faults. The results obtained are directly applied 
to the real artworks in this work with the assumption that the test sample contains a wide range of faults 
from very stronger ones to very weaker ones. Nevertheless, the measure can be improved by rigorous 
training on simulated faults of different shapes, sizes, depth and orientations which goes beyond the scope 
of current work. 

The second important parameter in the proposed algorithm is the number of levels till which the original 
space should be decomposed with the selected basis function. Wavelet decomposition level, when viewed 
in terms of filter banks, determines the frequency resolution of the proposed method and will thus depend 
on the nature of the defects that need to identified. Although, multiple factors may influence this level 
selection including the geometry and physical properties of the material, the two fundamental parameters 
that are considered for level selection are the depths and the diameters of the defects. We exploit the 


min 
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laboratory developed test sample (Fig. |2l) to establish a criterion for level selection by computing mutual 
information based similarity index between the original image and its subsequent low pass versions at 
different levels. 

The wavelet decomposition yields higher frequency information of the analyzed content in lower 
decomposition levels and lower frequency information in the higher decomposition levels. The faults that 
are closest to the surface and widest manifest themselves with the highest intensities in the acquired 
thermal images vis-a-vis other faults of smaller dimensions that are farther from the surface. As a result, 
stronger the fault in terms of its diameter and depth, the more likely it is to appear in the lower frequency 
bands of the spectrum and thus in the lower decomposition levels. Based on this assumption, while the 
smaller and weaker faults may be retrieved in some latest levels of decomposition, the stronger faults 
are likely to be associated with the background in the subsequent levels. The wavelet decomposition 
essentially involves successive application of two-channel (low pass and high pass) filtering of the 
approximation (low pass content) at each level with perfect reconstruction capability. Hence, when at a 
certain level L, information about the strongest fault has been completely detached from the approximation 
at that level, Al, the remaining fault information essentially goes into the detail Dl at that level. Any 
further decomposition performed after this level L would thus yield no additional information regarding 
the faults. We utilize a measure based on the mutual information for automating the level selection. 

Mutual information between two random variables allows measuring the similarity between them 
and the mutual information based measures have been extensively used as similarity index in image 
processing domain |[22]| . Il23l . f2A\ . |[25l . Ehl . ETl . A slight variant of the classical mutual information 
is the Regional Mutual Information (RMI) ll^ . RMI exploits localized regional information rather than 
the global information and has been demonstrated to be a more robust method for measuring the mutual 
information between two images. A higher RMI indicates a higher degree of dependence between the 
analyzed images. 

We compute the RMI between the raw thermal image, Yt and a subspace spanning the approximation at 
level L for every pulse in the PRBS sequence. The resulting RMI is normalized by the mutual information 
of the raw image and then temporally averaged to obtain Mavg{L) for a given basis function. The resulting 
mutual information (Fig. |7]l follows a monotonous trend as expected with reducing similarity index at 
each subsequent level. Mavg{^) is the highest value obtained at the first level of decomposition and the 
difference between the successive level values remains constant until level 6. However, from level 6 
onwards, this difference becomes relatively small and stays that way for subsequent levels. It can thus 
be deduced from this graph that level 6 represents the maximum decomposition level for the test sample 
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being analyzed. At this level, the most dominant fault has been completely detached from the background 
and the pictorial layers and hence the mutual information does not change significantly for higher levels 
of decomposition. This observation was ascertained by visual inspection of the approximation at level 6, 
where the Hole 1 was found to have been removed. An other way to argue this is to use the scale(length 
in X and y direction) of the analyzing wavelet functions. This can be estimated in the x direction by 
2^-^Nx and in the y-direction by the 2^~^Ny where L is the level of wavelet decomposition.We can 
see that any value of L beyond 6 would result in a very small length and can hence be thought of as 
noise and the L = lor2 result in big lengths and can be interpreted as capturing the background and 
illumination effects.As discussed in the Section III-AI the distance between the camera and the mural is 
fixed so the correspondence between the pixel in image and the real distance on the mural surface can 
be established via simple trigonometry. Hence a similar argument can be established to match different 
wavelet decomposition levels to size of faults and irregularities represented therein. The proposed level 
selection measure was tested for different basis functions excitation times (Tg) but similar trend (difference 
between successive levels) was observed albeit with minor variations in the parameter values as shown 
in the Figure |7] . The proposed measure presents the best compromise for the problem at hand and while 
it was developed based on laboratory generated defects, the obtained parameter (level of decomposition) 
can be utilized for real artworks having similar fault dimensions. At this point, it remains difficult to give 
a generic parameter, which caters for all types and dimensions of defect geometries. In a future work, 
it would be interesting to study the problem in the context of test sample characterization through its 
physical properties and to correlate with the findings of fhe proposed method. 

The framework discussed above allows an efficient selection of model parameters for the proposed 
algorithm. The useful subspace is then reconstructed by considering intermediate level details. The higher 
frequency details account for the noise, whereas the Lth level approximation caters for the background 
affects. In the next section, the results obtained by application of the method to the test sample and the 
real artwork are presented along with elaborate discussion on the revealed fault maps. 

IV. Results and Discussion 

A. Application to Test Sample 

Fig. |8(a)| illustrates one of the raw images, captured on the test sample (Fig. O, with the best visibility 
of smaller faults. It can be observed from this image that there is a strong influence of the illumination 
pattern on the acquired image. This pattern results from non-homogeneous illumination of the test sample 
by the two halogen lamps. An important processing task will thus be to suppress this pattern to reveal 
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Fig. 1. The normalized region mutual information between the raw data Yt and the approximation at level L. The mutual 
information monotonically decreases with increasing L before being stagnant after reaching L — 6. The rate of change becomes 
insignificant after L — 6, suggesting the complete removal of separable faults from the approximation at level 6. It can be 
observed that the trend is the same across all the excitation times (Te) 


TABLE III 

Cost function values (Eq.O for a set of wavelet basis function for different excitation times (Te), d?, 

TESTED TO DETERMINE THE OPTIMAL BASIS.THE OPTIMAL VALUES FOR DIFFERENT EXCITATION TIMES ARE IN BOLDFACE. 


Basis 

d&4 

dblO 

coif! 

coif 4: 

coif 5 

sym4 

symlQ 

bior2.2 

bior2.4 

biorS.S 

bior3.7 

bior3.9 

bior6.8 

rbio3.7 

rbioQ.8 

Te = 0.5s 

9.837 

9.583 

8.340 

9.346 

9.563 

8.894 

9.088 

8.658 

8.5910 

11.0385 

10.0823 

9.8274 

8.9697 

8.9650 

8.9288 

Te = Is 

8.9912 

8.9888 

7.6042 

8.383 

8.7456 

7.5807 

8.1286 

7.9584 

7.7010 

10.8070 

10.0097 

9.7387 

7.9985 

8.9699 

7.9613 

Te = 2s 

7.8232 

7.3634 

6.8023 

7.8996 

7.7596 

6.097 

7.2485 

6.922 

6.7302 

8.4880 

8.1136 

8.1599 

7.1465 

8.5845 

7.1146 

Te = 5s 

8.087 

7.951 

8.957 

8.239 

8.509 

7.761 

7.855 

8.403 

8.009 

8.903 

8.7 

8.683 

7.766 

9.1 

7.685 

Te = 10s 

8.9017 

8.9301 

7.8784 

8.4588 

8.6935 

7.9610 

8.2871 

7.9732 

7.9401 

10.5622 

9.6351 

9.3328 

8.1533 

8.8919 

8.1085 


more relevant information regarding the defects. The most dominant hole, Hole l(4mm depth, 10 mm 
diameter), is easily ohservahle on this raw image, while Hole 2 also appears on the raw data alheit with 
reduced intensity. However, all the other holes have been masked hy the strong illumination effect and 
the general texture of the test material. 

In the first instance, we consider application of the proposed method to the laboratory developed test 
sample (Fig. II. The first parameter that needs to be selected is the set of basis functions for data 
decomposition and reconstruction. In this study, we tested different standard wavelet basis, a subset of 
which is listed in Table [Till Adopting the criterion developed in Eq. |2j the proposed measure for wavelet 
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basis selection is also presented in Table |III] Although the analysis was performed for different excitation 
times, results and the fault maps for only one (Tg = 5s) are displayed (Figure [8]l due to the shortage of 
space. The following discussion hence discusses this one case. The cost function has comparable values 
for a number of different basis functions and the minimum value is achieved for the basis function 15, 
as shown in Table |III1 i.e., the reverse biorthogonal rbio 6.8. It may be argued that the selected basis 
is not uniquely the optimal basis, which will provide the best results in all the scenarios, but it has 
been derived systematically rather than through trial and error. All the subsequent analysis on both the 
laboratory developed sample and the real art works will be carried out using this optimal basis set. 

The application of the proposed algorithm to the test sample of Fig. [2] leads to the detection results 
illustrated in Fig. |8(b)[ The data are decomposed till level 6 and the useful subspace is reconstructed 
from the details obtained at level 3 to level 6. In order to highlight the efficacy of the proposed algorithm, 
a comparison with the previously developed SVD based approach ifTTl is performed. The SVD based 
subspace decomposition approach works on temporal signatures for each pixel and reconstructs the useful 
subspace after subtraction of a reference subspace, corresponding to the first singular value, and the high 
frequency noise subspace. The higher order statistics of skewness and kurtosis are then computed over 
the useful subspace temporally to obtain the final detection images demonstrated in Fig. |8(c)|8(d) The 
results of the proposed approach clearly outperform the results of the SVD based approach where none 
of the defects, except a few traces of Flole 1, appear in the detection results. A comparison with the raw 


data (Fig. |8(a)[ ) suggests that the SVD based approach was unsuccessful in achieving a good separation 
between the defects and the background. This can primarily be attributed to the inherent nature of 
SVD which is based on diagonalization using second order statistics. Moreover, the SVD based approach 
involves temporal processing of the data while considering spatial independence. The temporal evolutions 
of faulty and non-faulty pixels in the acquired image follow a very similar trend and while, there are 
subtle variations in their respective temporal behaviours, the difference is not significant to establish a 
pronounced distinguishability. The suppression of the subspace corresponding to the first singular value 
therefore suppress useful information related to the defects. 


A comparison of the results of the proposed approach (Fig. |8(b)[ ) with the raw dataset (Fig. |8(a)[ ) 
clearly shows that the proposed approach allows identification of most of the defects in the test material 
with a better spatial localization. In particular, of the total 12 holes, Floles 1-10 are visible now after 
processing through the proposed system. The faults corresponding to Floles 1 and 2 are the most dominant 
followed by Holes 3, 5, 6, 9 and 10. The Holes 4, 7 and 8 appear with relatively lesser intensities as 
compared to other faults due to their dimensions. Finally, the weakest faults Holes 11 and 12 are not 
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detected by our proposed approach. It is pertinent to recall that the thermal radiography is a technique, 
which works well closer to the surface and these weaker holes are probably too small and too far from 
the surface. Their non-detectability is thus more to do with the physical limitation of the acquisition 
process than the proposed processing scheme. Moreover, for a given depth, the intensity of the detected 
faults enhances with increasing diameters and similarly for a given diameter, the intensity reduces as the 
faults move farther from the material surface. Although, at this stage, it is not the interest to accurately 
quantify the defects in terms of their dimensions but the relative intensities of different holes in the final 
detection parameter is a potential indicator of the fault dimension (depth and diameter). This opens up 
the possibility of exploring indirect fault quantification oriented approach in the future. 

In order to give a qualitative measure for the utility of the proposed approach, we compare the Signal 
to Noise Ratio (SNR) for different defects in the raw data against the SNR in the processed detection 
results. In order to compute the SNR after the application of the proposed algorithm, the signal power for 
a given hole is computed in a window centered around that hole, whereas the noise power is computed 
by averaging the powers in different background windows not containing the faults, being the same with 
Bavg used in Eq. |2] This process can be represented by the equation (1^ . 

The value jopt is the optimal basis for the wavelet decomposition as described in the section njl-R2l 
and k is the hole number thus k = 1,2, ■ ■ ■ ,12. The SNR for the raw data (SNRraw) is calculated in a 
similar way by replacing the Ydet(win)(k,j-opt) with y{win)k as represented in (l3bl) . To formulate a 
measure of the performance of the proposed algorithm denoted SNRimprv we calculate the improvement 
in the SNR resulting from the application of the algorithm using the equation (1^ . 


SNR,,d = 


Ydet(wm 


(kj-opt) 


Bavg,jopt 

3^(win)(k) 


SNRraii, = 


B 


avg,jopt 


SNR, 


imprv 


= SNRy,d - SNRr 


(3a) 

(3b) 

(3c) 


These SNRimprv values for various excitation times are shown in Table HIl 

It can be observed that the SNRs for most of faults were very low in the raw data and the Hole 1, 
representing the most dominant hole, was almost at the same level as the background. The other holes 
scored progressively lower with respect to SNR in the raw data. A marked improvement is obtained 
after processing the data through the proposed algorithm as seen in the last column of Table |II1 Hole 
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1 now has an SNR of 9.7 dB whereas most of the other holes which were not visible in the raw data 
have been significantly enhanced to higher level SNRs. This clearly demonstrates the strength of the 
proposed approach to extract the useful information from the raw thermal images. Some of the lower 
holes, however, demonstrate lower SNRs but still show some visible traces in the detected image. The 
improvement in the SNR achieved by the application of the proposed methodology on the rest of the 
excitation times (Tg) has been included in the Table ini It demonstrates that the results obtained are quite 
similar to the case discussed above and that the proposed algorithm works effectively. 

B. Application on the real artwork 

After validation of the proposed method on the laboratory generated test sample, we focus our attention 
to the real artwork. All the parameters of the proposed algorithm including the wavelet basis and the 
decomposition levels for Mural 1 are kept same as those for the laboratory generated test sample. Fig. 
|9(a)| shows the best raw image available for Mural 1. The detection results obtained for the Mural 1 are 
shown in Fig. |9(b)| which demonstrate a very good detection of all faults. Faults C and D are highlighted 
with same intensities. This is very interesting since they are both at the same, 3 mm, depth. Fault A is 
highlighted in the bottom with a higher intensity. It is interesting to note that this fault manifests variable 
detection intensities in its proximity, a fact attributable to the varying depth of this defect from 3 mm to 
10 mm. Fault B, which is deeper at 5 mm, shows a lighter intensity which correspond to the center of 
the fault A. A region around the eyes in Mural 1 is also detected in the result, which is a false alarm in 
this case. The defect E is not detectable by the proposed approach, owing partially to its deep location 
(10 mm) from the surface of the material. However, a slight heterogeneity on the upper right part of the 
fault D might be attributable to this fault. Compared with the state-of-the art method (results shown in 
Figs. |9(c)| and |9(d)| ), the shapes of the faults are better identified, especially for the ones closest to the 
front side (C, D and B). Moreover, this systematic identification of faults is available through a single 
output of the proposed method. This is in contrast to the skewness, which identifies only the faults C 
and D and to the kurtosis, which identifies only the faults A and B, albeit with a poor shape. 

The proposed method was tested on other data sets as well with success. This non-invasive method, 
thus allows us to characterize the faults in artworks using a simple setup with a heating distributed over 
the time by a PRBS excitation and a strong processing strategy. The approach can be adapted to other 
types of materials with minimal modifications. 

The proposed algorithm is next applied to Mural 2 to obtain the detection results as shown in Fig. IIV-BI 
While, the raw image depicted a strong anomaly in the top left corner (Fig. IIV-BI ). the detection results 
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(e) Fault map obtained by proposed algorithm 

Fig. 8. The detection results for the laboratory generated test sample of FigO: (a) raw data Yt showing only Floles 1 and 2; 
(b) estimated Yuet with the proposed scheme, significantly improving identification of defects over the raw data; (c)-(d) the 
detection results of the SVD based approach showing the skewness and kurtosis, which do not reveal much information about 
the defects.(e) the fault map is calculated for the proposed algorithm by ignoring the edge artifacts and thresholding 
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(c) Skewness of the SVD based data reconstruction 


(d) Kurtosis of the SVD based data reconstruction 


Fig. 9. Detection parameter for the Mural 1: (a) Raw thermal image at an acquisition instant providing best visibility of faults 
in raw data; (b) Detection parameter of proposed scheme improving on the raw results; (c)-(d) Results of SVD based approach 
showing skewness and kurtosis. The contours mark the edges of the faults A to E (Fig. [3}. 4 defects are revealed by the 
proposed approach in (b) with some anomalies corresponding to the eye region of Mural 1. 


show three additional anomalies. In order to obtain a closer inspection of the identified anomalies, the 
temporal evolution of the the reconstructed suhspace is plotted in Fig. [TT] for four different regions 
including two principal anomalies on the right side (top and bottom), a background location and the 
anomaly on the top left. It was observed that the anomaly on the top left exhibits a very high temporal 
correlation with the excitation sequence, whereas the temporal correlations for the other anomalies are 
relatively lower. The results are then post adjusted by setting a correlation based threshold on the detection 
results to remove the anomalies that are highly correlated with the excitation sequence. The results after 
applying this correction are demonstrated in Fig. |10(c)[ The three anomalies originally identified are 
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■ 302 



- 0.2 


50 100 150 200 250 

(a) The raw data for Mural 2 


50 100 150 200 250 
(b) The fault map for Mural 2. 


50 100 150 200 250 


(c) The corrected detection after remov¬ 
ing correlated data. 


Fig. 10. The detection results for Mural 2:(a)The raw data for the Mural2 , gold foil is visible (b) the detection parameter of 
the proposed algorithm depicting four anomalies; (c) the correlation corrected image allowing removal of the anomaly on top 
left of (a) that was associated with a gold foil. 


retained in the detection parameter whereas the fourth anomaly has been subdued to a large extent. This 
fourth anomaly was actually due to a gold foil in that spatial vicinity which resulted in a high reflection 
and thus a higher correlation with the excitation sequence. The analysis thus allows us not only to detect 
the defects in the real artworks with a better separation from the background but also to identify the 
anomalies that may result from highly reflective surfaces. Out of these regions the one on the top left 
is an anomaly caused by a gold foil coating on the surface of the specimen being analysed. On closer 
inspection we find that the first point is not a conventional fault but is a piece of gold foil on the pictorial 
layer. This leads to its distinct behavior when we plot the time series specified by this point in Yrecon as 
shown in Figure [TT] The point visually shows marked similarity between the PRBS excitation sequence 
and the time-series Yrecon{W, 100). When we take the correlation of all the four time-series with the 
PRBS excitation sequence the result is depicted in Figure |1 1(b) Two points which yield a very high 
correlation are the T).econ(19,100) and the Yrecon{^^7, 140) the former being the gold foil and later 
the background.Based on the marked difference in the correlation of the different time-series we apply 
thersh-holding on the maximum value of the correlation coefficient to eliminate the effect of reflective 
materials (gold foil)and thus obtain an improved estimate of the fault location as shown in Figure |10(c)[ 


August 26, 2015 


DRAFT 





















24 



(a) Comparison of the Timeseries of Fault, Background 
and the Gold Foil obtained after proposed algorithm 


6 



200 250 300 350 400 450 500 

Time Samples 


(b) PRBS excitation sequence used for the Mural 2 
Analysis 

Fig. 11. Temporal evolution of the reconstructed subspace at four different spatial locations for the Mural 2 including three 
anomalies (top left, top right, top left) and a background location. The time series corresponding to the anomaly on top left 
exhibits a high correlation with the excitation sequence and this anomaly thus represents a highly reflective surface. 


V. Conclusion 

Non-destructive testing of materials is a highly developed field where different techniques are employed 
for automated characterization of test materials. One particular type of material that is highly sensitive to 
rupture is ancient artwork. Photo-radiometry can he employed to test these artworks, where the material 
is heated and the resulting thermal response is captured hy a thermal infrared camera. In this paper, a 
pseudo-radom binary sequence based excitation, as opposed to the classical pulse excitation, is employed 
which prevents excessive heating of the material thus reducing the risk of irreversible damage to the 
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artworks. Using simplified thermal models, the faults in the material should induce an irregularity in 
the acquired thermal image at their spatial locations. However, owing to different factors such as the 
background pictorial layer, illumination conditions, the measurement noise and the diffusion, the faults are 
not readily identifiahle on the raw acquired images. While, the faults of larger dimensions are visible in the 
raw data, they are at similar intensity levels as the background. The weaker faults are masked by the afore 
mentioned factors and are not visible at all in the raw data. Considering the spatial non-stationarity in the 
acquired raw data, we developed an algorithm based on the wavelet subspace decomposition. Appropriate 
selection criteria for selection of different parameters like the basis selection and the decomposition level 
selection were also proposed in this paper. The application on a laboratory generated test sample allowed 
establishing these parameters and the final subspace obtained represented the fault maps in the spatial 
domain, thereby allowing characterization of different faults. A significant enhancement in the SNR of 
the faults was obtained in the detected subspace as compared to the raw data. The method was also 
demonstrated to successfully identify faults in real artworks. It was also observed that the intensities of 
the detected fault maps are proportional to the fault dimensions (diameter and depth), a potential indicator 
for quantification of faults in future work. In brief, this work integrates a simple non-invasive material 
testing setup with a post-processing strategy, enabling an efficient spatial localization of faults in the 
material. 
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